
clear all

********************************************
* Estimate model for Table D1
********************************************	

local cell = "A1"
local append_replace = "replace"
local iter = 0

foreach var in 1 3 7{
local ++iter	
	
if `iter' == 2{
	local cell = "F1"
	local append_replace = "modify"
}

if `iter' == 3{
	local cell = "K1"
}

use "$path/Data/working.dta", clear
keep if sample == 1

**** Extract treatment effects
qui sum d1_`var' if treatment == 6
	local Y_C = r(mean)

* Overall Support (Table 2 Column 2)
reg d1_`var' ib6.treatment i.strata i.wave phone_survey survey_date bd1_`var', cluster(ent_id)
	local Y_LG = `Y_C' + _b[5.treatment]
	local Y_G = `Y_C' + _b[3.treatment]
	local Y_I = `Y_C' + _b[4.treatment]

* Treatment impacts on awareness of aid sharing OR not resenting refugees as aid recipients
gen aware_or_noresent = d3_3 == 1 | (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != . & d3_3 != .
reg aware_or_noresent ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]

*** Solve model
	
* Compute transformed support vars
foreach t in LG G I C{
local S_`t' = log(-log(1-`Y_`t''))
}

* Define LHS values
local Diff_S_LG = `S_LG' - `S_C'
local Diff_S_G = `S_G' - `S_C'
local Diff_S_I = `S_I' - `S_C'

* Create matrices
matrix Y = (`Diff_S_LG', `Diff_S_G', `Diff_S_I')'
matrix X = ( ///
    `Diff_A_LG', 1 , 1 \ ///
    `Diff_A_G', 1 , 0 \ ///
    `Diff_A_I', 0 , 0 ///
)

* Solve for coefficients
matrix b = syminv(X' * X) * X' * Y

scalar A = b[1,1]
scalar W = b[2,1]
scalar alpha_LG = b[3,1]

display "A = " A
display "W = " W
display "alpha_LG = " alpha_LG

*** Verify the solution
local Y_verify_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + W + alpha_LG))
disp "Y_verify_LG = "`Y_verify_LG'
disp "Y_LG = "`Y_LG'

local Y_verify_G = 1 - exp(-exp(`S_C' + A * `Diff_A_G' + W))
disp "Y_verify_G = "`Y_verify_G'
disp "Y_G = "`Y_G'

local Y_verify_I = 1 - exp(-exp(`S_C' + A * `Diff_A_I'))
disp "Y_verify_I = "`Y_verify_I'
disp "Y_I = "`Y_I'

* Benchmark 
local Diff_Y_LG = `Y_LG' - `Y_C'
local Diff_Y_G = `Y_G' - `Y_C'
local Diff_Y_I = `Y_I' - `Y_C'

*** Compute counterfactuals
local Diff_Y_nowealth_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + alpha_LG)) - `Y_C'
disp "Diff_Y_nowealth_LG = "`Diff_Y_nowealth_LG'

local Diff_Y_wealth_G = 1 - exp(-exp(`S_C' + W)) - `Y_C'
disp "Diff_Y_wealth_G = "`Diff_Y_wealth_G'

*** Output results
matrix app_d_output = (A, W, alpha_LG, `Diff_Y_nowealth_LG', `Diff_Y_wealth_G')
matrix colnames app_d_output = A_`var' W_`var' alpha_LG_`var' Diff_Y_nowealth_LG_`var' Diff_Y_wealth_G_`var' 

cd "$path/Data"
putexcel set app_d_output.xls, `append_replace'
putexcel `cell' = matrix(app_d_output), colnames
}

********************************************
* Bootstrap standard errors
********************************************
local constrain_count = 0
local cell ""
use "$path/Data/working.dta", clear
keep if sample == 1
set seed 123456789

forvalues i=1/2000 {
local row = `i' + 2
local cell = "A"+string(`row')

preserve
bsample, strat(treatment strata wave)

* Treatment impacts on awareness of aid sharing (Table 3 Column 3)
gen aware_or_noresent = d3_3 == 1 | (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != . & d3_3 != .
qui reg aware_or_noresent ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]

* Support Outcomes
foreach var in 1 3 7{

if `var' == 3{
local cell = "F"+string(`row')
}

if `var' == 7{
local cell = "K"+string(`row')
}

**** Extract treatment effects
qui sum d1_`var' if treatment == 6
	local Y_C = r(mean)

* Support (Table 2)
qui reg d1_`var' ib6.treatment i.strata i.wave phone_survey survey_date bd1_`var', cluster(ent_id)
	local Y_LG = `Y_C' + _b[5.treatment]
	local Y_G = `Y_C' + _b[3.treatment]
	local Y_I = `Y_C' + _b[4.treatment]

*** Solve model
	
* Compute transformed support vars
foreach t in LG G I C{
local S_`t' = log(-log(1-`Y_`t''))
}

* Define LHS values
local Diff_S_LG = `S_LG' - `S_C'
local Diff_S_G = `S_G' - `S_C'
local Diff_S_I = `S_I' - `S_C'

* Create matrices
matrix Y = (`Diff_S_LG', `Diff_S_G', `Diff_S_I')'
matrix X = ( ///
    `Diff_A_LG', 1 , 1 \ ///
    `Diff_A_G', 1 , 0 \ ///
    `Diff_A_I', 0 , 0 ///
)

* Solve for coefficients
matrix b = syminv(X' * X) * X' * Y

scalar A = b[1,1]
scalar W = b[2,1]
scalar alpha_LG = b[3,1]


if (scalar(A) < 0) {
	scalar A = 0 // constrain gamma_A to be non-negative
	scalar W = `Diff_S_G' // if gamma_A = 0, gamma_W = S_G - S_C
	scalar alpha_LG = `S_LG' - `S_G' // if gamma_A = 0, gamma_W = S_G - S_C and alpha_LG = S_LG - S_G
	local constrain_count = `constrain_count' + 1
}


* Benchmark 
local Diff_Y_LG = `Y_LG' - `Y_C'
local Diff_Y_G = `Y_G' - `Y_C'
local Diff_Y_I = `Y_I' - `Y_C'

*** Compute counterfactuals
local Diff_Y_nowealth_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + alpha_LG)) - `Y_C'
local Diff_Y_wealth_G = 1 - exp(-exp(`S_C' + W)) - `Y_C'

matrix app_d_output = (A, W, alpha_LG, `Diff_Y_nowealth_LG', `Diff_Y_wealth_G')
	
cd "$path/Data"
putexcel set app_d_output.xls, modify
putexcel `cell' = matrix(app_d_output)

}
restore
}


**gen sum stats from bootstrap results
cd "$path/Data"
import excel using app_d_output.xls, clear first

gen id = _n
reshape long A_ W_ alpha_LG_ Diff_Y_nowealth_LG_ Diff_Y_wealth_G_, i(id) j(var)
sort var id

foreach var in 1 3 7{
// bring in incorrect SEs (these will later be replaced through substitute)
eststo td1a_`var': mean A_ W_ alpha_LG_ if var == `var' & id != 1 
eststo td1b_`var': mean Diff_Y_nowealth_LG_ Diff_Y_wealth_G_ if var == `var' & id != 1

// store empirical confidence intervals and incorrect SEs in locals
local it = 0
foreach out in A_ W_ alpha_LG_ Diff_Y_nowealth_LG_ Diff_Y_wealth_G_{
local ++it
qui sum `out' if var == `var' & id != 1, d
local low = round(r(p5),0.01)
local low_fmt : display %4.2f `low'
local high = round(r(p95),0.01)
local high_fmt : display %4.2f `high'
local ci_`var'_`it' = "`low_fmt'"+","+"`high_fmt'"

qui reg `out' if var == `var' & id != 1
local se_`var'_`it' : display %6.5f _se[_cons]

}

// Recover parameter estimates and p-values from bootstrap results
local iter = 0
foreach outcome in A_ W_ alpha_LG_ Diff_Y_nowealth_LG_ Diff_Y_wealth_G_ {
	local++ iter
	qui summarize `outcome' if var == `var' & id == 1
		local par_est = r(mean)
	qui summarize `outcome' if var == `var' & id != 1
		local boot_mean = r(mean)
	matrix est_d`var'_`iter' = `par_est'
	qui count if var == `var' & id != 1 & `outcome' <= 0
	local share_low = r(N)/2000
	qui count if var == `var' & id != 1 & `outcome' >= 0
	local share_high = r(N)/2000
	matrix p_boot_d`var'_`iter' = 2*min(`share_low',`share_high')
}

matrix est_boot = (est_d`var'_1, est_d`var'_2, est_d`var'_3)
mat colnames est_boot = A_ W_ alpha_LG_
estadd matrix est_boot : td1a_`var'
matrix drop est_boot est_d`var'_1 est_d`var'_2 est_d`var'_3

matrix est_boot = (est_d`var'_4, est_d`var'_5)
mat colnames est_boot = Diff_Y_nowealth_LG_ Diff_Y_wealth_G_
estadd matrix est_boot : td1b_`var'
matrix drop est_boot est_d`var'_4 est_d`var'_5

matrix p_boot = (p_boot_d`var'_1, p_boot_d`var'_2, p_boot_d`var'_3)
mat colnames p_boot = A_ W_ alpha_LG_
estadd matrix p_boot : td1a_`var'
matrix drop p_boot p_boot_d`var'_1 p_boot_d`var'_2 p_boot_d`var'_3

matrix p_boot = (p_boot_d`var'_4, p_boot_d`var'_5)
mat colnames p_boot = Diff_Y_nowealth_LG_ Diff_Y_wealth_G_
estadd matrix p_boot : td1b_`var'
matrix drop p_boot p_boot_d`var'_4 p_boot_d`var'_5
}


label var A_ "Awareness Coefficient ($\gamma_A$)"
label var W_ "Wealth Term ($\gamma_W \Delta_W$)"
label var alpha_LG_ "Labeled Grant Fixed Effect ($\alpha_{LG}$)"
label var Diff_Y_nowealth_LG_ "Labeled Grant: No Wealth Effect"
label var Diff_Y_wealth_G_ "Grant: Wealth Effect"

** Significance stars need to be added to this table manually.

esttab td1a_1 td1a_3 td1a_7 using "$path/Output/Appendix_D/model.tex", label nostar nomtitles replace collabels(none) nolines nonumber noobs cells(est_boot(fmt(%9.2f)) se(par fmt(%6.5f)) p_boot(par([ ] ) fmt(%9.2f))) ///
substitute(\_ _ \$ $ "(`se_1_1')" "(`ci_1_1')" "(`se_1_2')" "(`ci_1_2')" "(`se_1_3')" "(`ci_1_3')" "(`se_3_1')" "(`ci_3_1')" "(`se_3_2')" "(`ci_3_2')" "(`se_3_3')" "(`ci_3_3')" "(`se_7_1')" "(`ci_7_1')" "(`se_7_2')" "(`ci_7_2')" "(`se_7_3')" "(`ci_7_3')") ///
prehead("\begin{table}[h!]	\centering	\footnotesize	\caption{Parameter Estimates} \label{tab:model} \begin{tabular}{l*{3}{>{\centering\arraybackslash}p{2.4cm}}} \toprule \toprule  & \shortstack{Supports\\Refugee\\Hosting} & \shortstack{Supports\\More\\Refugees} & \shortstack{Supports\\Right\\to Work} \\") ///
posthead("\cmidrule{2-4} \textit{Parameter Estimates} &&& \\ ") ///
prefoot("") ///
postfoot("")

esttab td1b_1 td1b_3 td1b_7 using "$path/Output/Appendix_D/model.tex", label nomtitles append collabels(none) nolines nonumber noobs cells(est_boot(fmt(%9.2f)) se(par fmt(%6.5f)) p_boot(par([ ] ) fmt(%9.2f))) ///
substitute(\_ _ \$ $ "(`se_1_4')" "(`ci_1_4')" "(`se_1_5')" "(`ci_1_5')" "(`se_3_4')" "(`ci_3_4')" "(`se_3_5')" "(`ci_3_5')" "(`se_7_4')" "(`ci_7_4')" "(`se_7_5')" "(`ci_7_5')") ///
prehead("") ///
posthead("\textit{Counterfactual Treatment Effects} &&& \\ ") ///
prefoot("") ///
postfoot("\bottomrule \bottomrule \multicolumn{4}{p{0.85\linewidth}}{\footnotesize See Appendix \ref{sec:structural} for estimation details. $ \gamma_A$, $ \gamma_W \Delta_W$, and $ \alpha_{LG}$ are obtained by solving Equation \ref{eqn:identification} under the constraint $ \gamma_A > 0$. \textit{Labeled Grant: No Wealth Effect} shows the counterfactual treatment effect of Labeled Grant when $ \gamma_W = 0$ estimated using Equation \ref{eqn:support}. \textit{Grant: Wealth Effect} shows the counterfactual treatment effect of Grant Only when $ \gamma_A = 0$ estimated using Equation \ref{eqn:support}. Bootstrap 90\% confidence intervals in parentheses show the 5th and 95th percentiles from 2,000 simulations. $ p $-values in brackets test the two-sided hypothesis that the coefficient or treatment effect equals zero, and are estimated by doubling the largest $\alpha$ for which the bootstrap $1-\alpha$ confidence interval includes zero. \sym{*} \(p<0.1\), \sym{**} \(p<0.05\), \sym{***} \(p<0.01\).}  \end{tabular} \\ \end{table}%")


estimates drop td1a_1 td1a_3 td1a_7 td1b_1 td1b_3 td1b_7

cd "$path/Code"
erase "$path/Data/app_d_output.xls"



*------------------------------------------------------------------------------------
* Footnote 48
*------------------------------------------------------------------------------------

disp "`constrain_count' estimates out of 6000 were constrained to gamma_A = 0"


*------------------------------------------------------------------------------------
* Footnote 44: Model solution when estimating awareness and resentment separately.
*------------------------------------------------------------------------------------


clear all
	
use "$path/Data/working.dta", clear
keep if sample == 1

**** Extract treatment effects
qui sum d1_1 if treatment == 6
	local Y_C = r(mean)

* Overall Support (Table 2 Column 2)
reg d1_1 ib6.treatment i.strata i.wave phone_survey survey_date bd1_1, cluster(ent_id)
	local Y_LG = `Y_C' + _b[5.treatment]
	local Y_G = `Y_C' + _b[3.treatment]
	local Y_I = `Y_C' + _b[4.treatment]

* Treatment impacts on awareness of aid sharing (Table 4 Column 3)
gen aware_or_noresent = d3_3 == 1 | (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != . & d3_3 != .
reg aware_or_noresent ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]
	
reg d3_3 ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]

gen noresent = (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != .
reg noresent ib6.treatment i.strata i.wave phone_survey survey_date, cluster(ent_id)
	local Diff_R_LG = _b[5.treatment]
	local Diff_R_G = _b[3.treatment]
	local Diff_R_I = _b[4.treatment]

*** Solve model
	
* Compute transformed support vars
foreach t in LG G I C{
local S_`t' = log(-log(1-`Y_`t''))
}

* Define LHS values
local Diff_S_LG = `S_LG' - `S_C'
local Diff_S_G = `S_G' - `S_C'
local Diff_S_I = `S_I' - `S_C'

* Create matrices
matrix Y = (`Diff_S_LG', `Diff_S_G', `Diff_S_I')'
matrix X = ( ///
    `Diff_A_LG', `Diff_R_LG' , 1 \ ///
    `Diff_A_G', `Diff_R_G' , 1 \ ///
    `Diff_A_I', `Diff_R_I' , 0 ///
)

* Solve for coefficients
matrix b = syminv(X' * X) * X' * Y

scalar A = b[1,1]
scalar R = b[2,1]
scalar W = b[3,1]

if (scalar(A) < 0) {
	scalar A = 0 // constrain gamma_A to be non-negative
	scalar W = `Diff_S_G' // if gamma_A = 0, gamma_W = S_G - S_C
	scalar alpha_LG = `S_LG' - `S_G' // if gamma_A = 0, gamma_W = S_G - S_C and alpha_LG = S_LG - S_G
}

display "A = " A
display "R = " R
display "W = " W

*** Verify the solution
local Y_verify_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + R * `Diff_R_LG' + W))
disp "Y_verify_LG = "`Y_verify_LG'
disp "Y_LG = "`Y_LG'

local Y_verify_G = 1 - exp(-exp(`S_C' + A * `Diff_A_G' + R * `Diff_R_G' + W))
disp "Y_verify_G = "`Y_verify_G'
disp "Y_G = "`Y_G'

local Y_verify_I = 1 - exp(-exp(`S_C' + A * `Diff_A_I' + R * `Diff_R_I'))
disp "Y_verify_I = "`Y_verify_I'
disp "Y_I = "`Y_I'

* Benchmark 
local Diff_Y_LG = `Y_LG' - `Y_C'
local Diff_Y_G = `Y_G' - `Y_C'
local Diff_Y_I = `Y_I' - `Y_C'

*** Compute counterfactuals
local Diff_Y_wealth_G = 1 - exp(-exp(`S_C' + W)) - `Y_C'
disp "Diff_Y_wealth_G = "`Diff_Y_wealth_G'

local Diff_Y_nowealth_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + R * `Diff_R_LG')) - `Y_C'
disp "Diff_Y_nowealth_LG = "`Diff_Y_nowealth_LG'

*** Output results
matrix app_d_output_alt = (A, R, W, `Diff_Y_wealth_G', `Diff_Y_nowealth_LG')
matrix colnames app_d_output_alt = A R W Diff_Y_wealth_G Diff_Y_nowealth_LG

cd "$path/Data"
putexcel set app_d_output_alt.xls, replace
putexcel A1 = matrix(app_d_output_alt), colnames

********************************************
* Bootstrap standard errors
********************************************
local cell ""
use "$path/Data/working.dta", clear
keep if sample == 1

forvalues i=1/2000 {
local row = `i' + 2
local cell = "A"+string(`row')

preserve
bsample, strat(treatment strata wave)

**** Extract treatment effects
qui sum d1_1 if treatment == 6
	local Y_C = r(mean)

* Overall Support (Table 1 Column 2)
qui reg d1_1 ib6.treatment i.strata i.wave phone_survey survey_date bd1_1, cluster(ent_id)
	local Y_LG = `Y_C' + _b[5.treatment]
	local Y_G = `Y_C' + _b[3.treatment]
	local Y_I = `Y_C' + _b[4.treatment]

* Treatment impacts on awareness of aid sharing (Table 3 Column 3)
gen aware_or_noresent = d3_3 == 1 | (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != . & d3_3 != .
qui reg aware_or_noresent ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]
	
qui reg d3_3 ib6.treatment i.strata i.wave phone_survey survey_date bd3_3, cluster(ent_id)
	local Diff_A_LG = _b[5.treatment]
	local Diff_A_G = _b[3.treatment]
	local Diff_A_I = _b[4.treatment]

gen noresent = (e2_aid7 >= 3 & e2_aid7 != .) if e2_aid7 != .
qui reg noresent ib6.treatment i.strata i.wave phone_survey survey_date, cluster(ent_id)
	local Diff_R_LG = _b[5.treatment]
	local Diff_R_G = _b[3.treatment]
	local Diff_R_I = _b[4.treatment]

*** Solve model
	
* Compute transformed support vars
foreach t in LG G I C{
local S_`t' = log(-log(1-`Y_`t''))
}

* Define LHS values
local Diff_S_LG = `S_LG' - `S_C'
local Diff_S_G = `S_G' - `S_C'
local Diff_S_I = `S_I' - `S_C'

* Create matrices
matrix Y = (`Diff_S_LG', `Diff_S_G', `Diff_S_I')'
matrix X = ( ///
    `Diff_A_LG', `Diff_R_LG' , 1 \ ///
    `Diff_A_G', `Diff_R_G' , 1 \ ///
    `Diff_A_I', `Diff_R_I' , 0 ///
)

* Solve for coefficients
matrix b = syminv(X' * X) * X' * Y

scalar A = b[1,1]
scalar R = b[2,1]
scalar W = b[3,1]

if (scalar(A) < 0) {
	scalar A = 0 // constrain gamma_A to be non-negative
	scalar W = `Diff_S_G' // if gamma_A = 0, gamma_W = S_G - S_C
	scalar alpha_LG = `S_LG' - `S_G' // if gamma_A = 0, gamma_W = S_G - S_C and alpha_LG = S_LG - S_G
}

* Benchmark 
local Diff_Y_LG = `Y_LG' - `Y_C'
local Diff_Y_G = `Y_G' - `Y_C'
local Diff_Y_I = `Y_I' - `Y_C'

*** Compute counterfactuals
local Diff_Y_wealth_G = 1 - exp(-exp(`S_C' + W)) - `Y_C'
local Diff_Y_nowealth_LG = 1 - exp(-exp(`S_C' + A * `Diff_A_LG' + R * `Diff_R_LG')) - `Y_C'

matrix app_d_output_alt = (A, R, W, `Diff_Y_wealth_G', `Diff_Y_nowealth_LG')
	
cd "$path/Data"
putexcel set app_d_output_alt.xls, modify
putexcel `cell' = matrix(app_d_output_alt)

restore
}

**gen sum stats from bootstrap results
cd "$path/Data"
import excel using app_d_output_alt.xls, clear first

sum Diff_Y_nowealth_LG if _n == 1
local par_est = r(mean)
sum Diff_Y_nowealth_LG if _n != 1
local boot_mean = r(mean)
count if _n > 1 & Diff_Y_nowealth_LG <= 0
	local share_low = r(N)/2000
count if _n > 1 & Diff_Y_nowealth_LG >= 0
	local share_high = r(N)/2000
local p_boot = 2*min(`share_low',`share_high')

* Footnote 44
disp "Effect of Labeled Grant Net of Wealth Effects (Alt Model) = `par_est'"
disp "p-value on that estimate = `p_boot'"

erase "$path/Data/app_d_output_alt.xls"
cd "$path/Code"

*------------------------------------------------------------------------------------
* Footnote 45: Simulate error of approximation in App D Eq 2
*------------------------------------------------------------------------------------

clear all

use "$path/Data/working.dta", clear
keep if sample == 1 & wave_midline == 1

local gamma_A = 1.1361524
local gamma_W = .10280185

qui sum bd1_1
local mu = log(-log(1-r(mean)))
gen epsilon = `mu' - 7.35 * log(-log(runiform())) // calibrate theta to match mean support

// Initialize a matrix to store results
matrix results = J(2000, 1, .)

// Perform the bootstrap sampling and computation
forvalues i = 1/2000 {
	
    // Draw a bootstrap sample
    preserve
    bsample
    
	gen support = epsilon > - `gamma_A' * bd3_3 - `gamma_W' * bd9_2
	qui sum bd3_3
		local aware_mean = r(mean)
	qui sum bd9_2
		local wealth_mean = r(mean)	
	gen support_proxy = epsilon > - `gamma_A' * `aware_mean' - `gamma_W' * `wealth_mean'
	qui sum support
		local support_mean = r(mean)
	qui sum support_proxy    
		local support_proxy_mean = r(mean)
    
    // Store the result in the matrix
    matrix results[`i', 1] = abs(`support_proxy_mean' - `support_mean') / `support_mean'
    
    // Restore the original dataset
    restore
}

// Compute the mean of the first column of the results matrix
mata : st_matrix("mean_error", colsum(st_matrix("results"))/2000)

local mean_error = mean_error[1,1]

*------------------------------------------------------------------------------------

// Footnote 44
disp "Effect of Labeled Grant Net of Wealth Effects (Alt Model) = `par_est'"
disp "p-value on that estimate = `p_boot'"

// Footnote 45: Display the mean
display "Mean absolute error = " `mean_error'

// Footnote 48
disp "`constrain_count' estimates out of 6000 were constrained to gamma_A = 0"
